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This paper describes an isobaric semi-grand canonical ensemble Monte Carlo scheme for the accu- 
rate study of phase behaviour in ternary fluid mixtures under the experimentally relevant conditions 

of prescribed pressure, temperature and overall composition. It is shown how to tune the relative 
chemical potentials of the individual components to target some requisite overall composition and 
how, in regions of phase coexistence, to extract accurate estimates for the compositions and phase 
fractions of individual coexisting phases. The method is illustrated by tracking a path through the 
composition space of a model ternary Lennard-Jones mixture. 



I. INTRODUCTION 

Fluids comprising multiple distinct components or 
'species' are pervasive in chemistry, physics and chemical 
engineering where they feature in contexts as diverse as 
chemical solutions, fuels, lubricants, microemulsions and 
cleaning products. Understanding the phase behaviour 
of such mixtures is important for exploiting their novel 
properties, for example, by ensuring a lubricant remains 
fluid under prescribed working conditions, or in designing 
processes for separation of hydrocarbons. 

In principle, the phase behaviour of a fluid mixture can 
be determined experimentally. However, the parameter 
space of composition, temperature and pressure is typi- 
cally large, rendering this costly in both time and money. 
Accordingly, there is much interest in efficiently predict- 
ing the phase behaviour of multicomponent fluids using 
computer simulation. In attempting to do so, it is rea- 
sonable to align oneself with the experimental situation 
in which (typically) one fixes all but one or two of the 
system parameters and varies the remaining ones within 
some bounds. Thus, for instance, in a ternary mixture 
one might set the temperature and pressure, and vary 
the composition of the mixture in a prescribed way, eg. 
by varying the relative concentration of a pair of com- 
ponents. Thus one explores the phase behaviour along a 
particular path in composition space. 

In this paper we shall consider a computational strat- 
egy for determining what happens if such a pathway in- 
tersects a region of phase coexistence. When this occurs 
the system will separate into two or more phases, and we 
seek to determine -accurately- their compositions and 
relative quantities accurately for all points on the path. 



II. PHASE SEPARATION AND 
FRACTIONATION 

The overall composition of a three component 
(ternary) mixture is specified in terms of the set of con- 
centrations of the individual components, .x^'^', a;2°'', X3°'* 
or, more succinctly, with i = 1 . . . 3. Here 



x\"> = NJN , (1) 

with Ni the total number of particles of species i in 
the system and X]?=i ^» ^^c total number of parti- 
cles. Clearly overall concentrations sum to unity, ie. 
'^i^i if ^ = 1, so that only two concentrations need be 
specified to define the composition. 

If for some prescribed overall composition, tempera- 
ture T, and pressure p, the system phase separates, the 
composition of the coexisting phases will, in general, dif- 
fer from the overall composition-a phenomenon termed 
'fractionation'. However conservation of particles implies 
that the concentrations of the individual phases are re- 
lated to the overall composition by a generalized lever 
rule. Let us assume that there are m coexisting phases, 
then 

m 
7=1 

Here x^^^ is the concentration of component i in phase 
7, defined as 

^l"^ = A^i"VE^i'^ (3) 

i=l 

while ^('•'^ is the "phase fraction" of phase 7 ie. the frac- 
tion of all particles that are found in phase 7: 

y-3 ^(7) 
^(7) ^ i (4) 

Of course phase fractions sum to unity, so ^'^'^^ = 1- 

It follows that in order to specify the coexistence prop- 
erties of the system at some {a;|'^'*}, one must determine 

{a;!''''} and ^^'^''^ for each phase present. Before describing 
how this determination can be made, we shall expand on 
the concept of a path in the space of overall composition. 
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III. PHASE DIAGRAM PATHWAYS 

Let us illustrate what we mean by a phase diagram 
pathway in the context of a ternary mixture. At fixed 
temperature and pressure a possible phase diagram for 
such a system is shown in Fig. 1 in the form of a Gibbs 
phase triangle [1] . The diagram includes a region of two 
phase coexistence with tie lines that shrink to zero at a 
critical point. These tie lines connect points that repre- 
sent the compositions of the coexisting phases. Any point 
of overall composition within the coexistence region will 
divide a tie line into two segments, the relative lengths of 
which yield the phase fractions of the coexisting phases. 
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FIG. 1. Schematic of an exemplary phase diagram of a ternary 
mixture at fixed temperature and pressure, represented in 
terms of a Gibbs phase triangle as described in the text. 

Also shown in Fig. 1 is an experimental pathway in 
composition space that passes through the coexistence 
region (dotted line). Existing simulation approaches are 
not well suited to determining coexistence properties 
along such a path. Instead they are tailored to the task 
of determining the boundary of the coexistence region it- 
self and the associated tie lines (see eg. [2-4]). Of course, 
from this information one can in principle infer the coex- 
istence properties (ie. the compositions and phase frac- 
tions) along the experimental path by graphical means- 
though the procedure may be of questionable accuracy. 
However, such a graphical solution is impracticable for 
mixtures of more than three components because com- 
pact representations of phase behaviour for such systems 
do not exist. Moreover, even for binary and ternary sys- 
tems, a straightforward graphical procedure is only fea- 
sible when both pressure and temperature are fixed: if 
the experimental path involves changes in temperature 
or pressure in addition to compositional changes, it be- 
comes necessary (with existing simulation methods) to 
map a portion of the full phase coexistence surface in 
order to infer the behaviour along the (one-dimensional) 
path of interest, ie. one needs to solve a higher dimen- 



sional problem than the one of actual interest. 

Thus it is of interest to develop new simulation strate- 
gies that are potentially more direct, accurate and flexible 
that those currently available and which may generalize 
to multicomponent mixtures and arbitrary paths through 
the phase diagram. In the current paper we take a step 
towards this goal by showing how one can track an arbi- 
trary path in the space of overall composition of a ternary 
mixture. 



IV. EXISTING SIMULATION APPROACHES 

Let us first briefly summarize principal existing simu- 
lation methodologies for obtaining the phase behaviour 
of ternary fluids. Two phase coexistence in a Lennard- 
Jones mixture was studied by Escobedo using Gibbs- 
Duhcm integration (GDI) in the semi-grand canonical 
ensemble [3]. GDI requires independent knowledge of 
a coexistence state point to bootstrap the integration, 
but once going it follows the coexistence binodal through 
composition space. Whilst this delivers more or less com- 
plete information on the phase behaviour (modulo inte- 
gration errors), it doesn't give direct information on the 

coexistence properties (ie. the {x^''^''} and ^(f)) along a 
particular experimentally relevant path in the phase dia- 
gram. Instead this behaviour has to be inferred from the 
knowledge of the binodal and the tie lines, as described 
in the previous section. 

Errington [4] has described both a semi-grand canoni- 
cal and a full grand canonical Monte Carlo (MC) scheme 
for studying fluid phase behaviour in multicomponent 
mixtures. This is a potentially powerful approach, as- 
pects of which we shall draw upon below. Nevertheless, 
it too is tailored for determining the coexistence binodal 
itself and without the extensions that we shall describe, 
is ill-suited to the task of tracking a particular experi- 
mental path and accurately determining the associated 
coexistence properties. 

In this respect Gibbs ensemble Monte Carlo (GEMC) 
is more flexible, since one can simply adjust the over- 
all densities of the components. However, as soon as 
the path to be followed passes close to the coexistence 
boundary (so that one of the phase fractions becomes 
small), the volume of one of the GEMC boxes also be- 
comes small and pronounced finite-size effects arc to be 
expected. The method can therefore only provide re- 
liable results for state points that are well within the 
coexistence region, leading to the same problems experi- 
enced by GDI and Errington's approach. The constant 
pressure GEMC technique has been applied by Tsang et 
al [2] to obtain the coexistence properties of the same 
ternary mixture that was studied by Escobedo [3]. 

Finally we mention a more bespoke technique for 
studying phase behaviour in mixtures, namely the pseudo 
ensemble method of Vrabec and Hasse[5]. This finds co- 
existence properties for prescribed choices of temperature 
and the composition of the liquid, which allows for cal- 
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culations of the binodal, but docs not seem capable of 
determining phase coexistence properties for prescribed 
values of the overall composition. 



V. METHOD 
A. Isobciric semi-grand canonical ensemble 

We turn now to the present method, which is in fact an 
adaptation of a fully grand canonical MC approach orig- 
inally developed for estimating cloud and shadow points 
in the phase diagrams of polydispcrsc fluids [6, 7]. Here, 
however, we shall work within the isobaric semi- grand 
canonical ensemble (SGCE). Within this framework the 
system volume V. the energy, and the instantaneous 
number densities of the species all fluctuate [8], while the 
particle number A'', pressure p, temperature T, and a set 
of chemical potentials {/ii} are all prescribed. The latter 
are measured with respect to an arbitrarily chosen refer- 
ence species, ie. for a ternary mixture there are actually 
only two independent chemical potential differences. 

Operationally, the sole difference between the iso- 
baric, semi-grand-canonical ensemble and the familiar 
constant-(A'',p, T) ensemble [9] is that one implements 
MC updates that select a particle at random and at- 
tempts to change its species label «, say, to another 
species i', chosen randomly. This proposal is accepted 
or rejected with a Metropolis probability controlled by 
the change in the internal energy and the chemical po- 
tential [10]. This probability reads: 

Pace min [1, exp (-[A$ + Hi - Hi']/T)] , 

where A$ is the internal energy change associated with 
the relabeling operation. 



B. Tracking a composition space pathway 

The SGCE is the appropriate ensemble for studying 
phase coexistence in multicomponent mixtures for two 
reasons. Firstly it realizes the common experimental sce- 
nario of flxed pressure and temperature; and secondly it 
provides for fluctuations in the densities of the individual 
species on the scale of the system size. This latter feature 
is crucial to allow for separation into differently fraction- 
ated phases. However, in order to track a composition 
space pathway, one first needs to bootstrap the tracking 
procedure by determining the set of chemical potentials 
that target some state point on the pathway. In what 
follows we describe how this can be achieved straightfor- 
wardly in two situations, namely (i) when the pathway 
passes through a one phase region of the phase diagram, 
or (ii) it intersects an axis of the ternary phase diagram, 
ie. coincides with the limit of a binary system. 



1. Bootstrapping the method in the one phase region 

The task is to determine, for some given T and p, the 
set of chemical potentials {/Uj} that yield some prescribed 

set of "target" concentrations {a;^°^}. By this we mean 
that the ensemble averages of the fluctuating overall con- 
centrations {xi} matches {a;-°^}. In a single phase re- 
gion, the non-equilibrium potential refinement (NEPR) 
scheme [11] enables the efficient iterative determination 
of {/Xi} from a single simulation run, and without the 
need for initial guesses. To achieve this, the {/Xj} are 
continually updated (in the course of a simulation run) 
in such a way as to minimize the deviation of the instan- 
taneous concentrations {xi} from the target. This pro- 
cedure realizes a non-equilibrium steady state for which 
{xi} = {xf^^}. However, since tuning model parameters 
'on-the-fly' in this manner violates detailed balance, the 
{/Uj} thus obtained is not the equilibrium solution one 
seeks. Nevertheless by performing a scries of iterations 
in which the degree of modification made to {/ii} at each 
step is successively reduced, one can drive the system 
towards equilibrium, thereby obtaining the correct {/ij}. 
Full details of the procedure are provided in [11]. An al- 
ternative scheme for determining chemical potentials has 
recently been described by Malasics et al [12]. 

Once {/ii} has been determined for a state point on the 
composition space path of interest, one can step along 
this path using histogram extrapolation techniques [13], 
without having to apply the NEPR technique again. The 
basic idea is measure fluctuating quantities at the initial 
state point on the path and reweight them to determine 
the {/ij} corresponding to a nearby state point on the 
path. The resiflting {/ii} are used to perform a fresh 
SGCE simulation at the new state point, and the process 
is repeated. 



2. Tracking the path into a two phase region 

Suppose now that the path being followed through 
composition space enters a two phase region. The bound- 
aries of this region can be located via extended sampling. 
The idea is to employ standard Multicanonical biasing 
techniques [14] to extend the range sampled by some 
order parameter. The precise choice of order parame- 
ter is not crucial, what matters is that it distinguishes 
clearly between the coexisting phases. In the case of a 
liquid-vapour transition, the system volume V is a suit- 
able order parameter and we monitor the order parame- 
ter probability distribution function (pdf), P{y). Close 
to the coexistence boundary this pdf exhibits a second 
much smaller peak representing contributions from the 
incipient phase. The appearance of the new phase means 
that fractionation will start to occur and in order to re- 
main on the desired path in overall composition space, 
one needs to account for this. To do so involves determin- 
ing values for the {/ij} that yield the prescribed overall 
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composition across the coexisting phases. At the same 
time, one wishes to determine accurately the {a;^''^} and 
notwithstanding the fact that the relative quantity 
of the incipient phase is close to zero. 

All this can be achieved via the following strategy. For 
given choices of , T, and p, one tunes {fii} and the 

^('•'^ iteratively within a histogram extrapolation frame- 
work, such as to simultaneously satisfy both a generalized 
lever rule and equality of the probabilities of occurrence 
of the phases, i.e. 



,(0)^^^(,)^(.)^ (5a) 

7=1 

£ = 0. (5b) 

In the first of these constraints. Eq. (5a), the ensem- 
ble averaged concentrations x'^'^^ are assigned by av- 
eraging only over configurations belonging to the re- 
spective phase[15]. The deviation of the weighted sum 
Xi = from the target x\^^ is conveniently 

quantified by a 'cost': 

A^^|S,-a.fM- (6) 
In the second constraint, Eq. (5b), 

7=1 ^ ^ 

provides a measure of the extent to which the probability 
of each phase occuring, tt^''^ , is equal for each of the m 
coexisting phases. Imposing this equality ensures that 
phase properties are measured under conditions of equal 
pressure and hence that finite-size errors in these proper- 
ties are exponentially small in the system volume [6, 16]. 
In practice the relative probabilities of the phases n^^^ is 
measured from the individual peak weights of the order 
parameter pdf. 

The determination of {/Uj} and ^^^^ such as to satisfy 
Eqs. (5a) and (5b) proceeds iteratively thus: 

1. Guess initial values of the S,^''^ corresponding to the 
chosen overall composition. Usually if starting near 
the boundary of the coexistence region, ^(^) for the 
incipient phase will be close to zero. 

2. Tune {fii} (within the histogram extrapolation 
scheme) such as to minimize A. 

3. Measure the corresponding value of £. 

4. if £ < tolerance, finish, otherwise vary ^(t) (within 
the histogram extrapolation scheme) and repeat 
from step 2. 



For the common case of two-phase coexistence, = 
1 — ^^^^ and the minimization of £ with respect to vari- 
ations in ^('•'^ is a one-parameter minimization which is 
easily automated using standard algorithms such as the 
"Brent" routine described in Numerical Recipes [17]. In 
step 2 the minimization of A with respect to variations 
in {ni} is most readily achieved [18] using the following 
simple iterative scheme for {fit}: 

/x^ = rt + aln(^j . (8) 

This update is applied simultaneously to each of the set of 
chemical potentials {/ii}, and thereafter the set is shifted 
so that /xi = 0, where species 1 is the chosen reference 
species. The quantity < a < 1 appearing in Eq. (8) 
is a damping factor, the value of which may be tuned to 
optimize the rate of convergence. In order to ensure that 
finite-size effects are exponentially small in the system 
size it is necessary to iterate to a high tolerance in step 
4; we chose this to be 10~^^. 

The values of ^^'^^ resulting from the; application of the 
above procedure are the desired phase fractions at the 
nominated {.x^'^'}. Properties of the individual phases 
such as their average number density, energy and concen- 
trations are obtainable by accumulating, for each phase 
separately, distributions of the respective quantity. 

Having satisfied the above procedure for the initial co- 
existence state point, one can track the path of interest 
deeper into the two phase region with the aid of his- 
togram extrapolation. One simply repeats the procedure 
but with target concentrations that correspond to a state 
point somewhat further along the path. Once the corre- 
sponding chemical potentials have been estimated, a new 
SGCE simulation will provide new data which in turn 
can be extrapolated yet further along the path. In this 
way one steps along the path obtaining the coexistence 
properties as one goes. 



3. Bootstrapping in the two phase region 

If the composition space pathway of interests termi- 
nates on an axis of the ternary phase diagram (ie. co- 
incides with the limit of a binary system) which fur- 
thermore, lies within a region of two phase coexistence, 
the tracking procedure can be bootstrapped straightfor- 
wardly at this state point. The strategy is first to elim- 
inate the third component by setting its chemical po- 
tential to a large negative value. Thereafter, tuning the 
difference in chemical potential between the remaining 
components allows one to scan along the coexistence tie 
line which itself necessarily coincides with the axis of 
the phase diagram (cf. Fig. 1). To achieve this one 
first links the phases via Multicanonical preweighting (see 
Refs. [19, 20] for strategies that achieve this), before ap- 
plying the methodology outlined above to determine the 
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chemical potentials and ^'•'^^ values corresponding to the 
terminal point on the pathway. Thereafter one increases 
the chemical potential of the third component so that it 
has a small but finite concentration, which then facili- 
tates tracking the path of interest via histogram extrap- 
olation, as outlined above. 



VI. ILLUSTRATIVE RESULTS 

We illustrate our procedure with results for a ternary 
Lennard Jones mixture, defined by 




/ \ 12 






(^) 







(9) 




where 



with 



o-22/o-ii = 0.75; 
e22/eii = 0.75; 



o-as/o-ii = 0.5 ; 
ess/en = 0.5 . 



FIG. 2. The thick black line indicates the path in overall 
(10) composition tracked in the simulations. Symbols denote the 
compositions of the coexisting phases (crosses are for the gas; 
filled circles are liquid), observed along this path. Statistical 
uncertainties do not exceed the symbol sizes. A selection of 
representative tie lines are also shown (dashed lines). At a;^' = 
0.409(1), the path leaves the region of phase coexistence. 



The potential was truncated at = S.OcTij and no tail 
correction was applied. 

We have studied the phase behaviour of this model 
for a system of = 256 particles at the reduced tem- 
perature T* = kBT/eu = 1.0 and the reduced pres- 
sure p* — pcTii/eii = 0.2. A similar system has been 
studied previously by Tsang et al [2], Escobedo[3] and 
Errington[4], although unlike the present work these 
studies utilized a correction for the potential truncation, 
so their results are not directly comparable with ours. 

The path we chose to study is one in which we fixed 
x^i ^ = 0.41 and varied X2 in the range [0,0.59]. The 
choice for 2:2'^'' then parameterizes the location along 
the path. The tracking procedure was bootstrapped at 
x^2 ^ = which corresponds to a binary system compris- 
ing components 1 and 3. As also observed by Tsang et 
al [2], vapor-liquid coexistence occurs in this limit. We 
therefore employed the bootstrapping approach outlined 
in sec. VB3 which yielded the value C^^) = 0.492(2) for 
the liquid phase fraction. Once bootstrapped, we per- 
formed a further run in which a small but finite value of 
exp(^2/r) was used in order to generate fluctuations in 

x^2 \ thus facilitating tracking of the path of interest in 
the manner described in Sec. VB2. 

Our results for the phase behaviour are summarized in 
fig. 2 which shows the path in overall composition (thick 
black line) as well as the compositions of the coexisting 
phases observed along this path (symbols). The corre- 
sponding measurements of the fraction of liquid, as 
a function of the path parameter x'^'^ are shown in Fig. 3. 
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FIG. 3. Measured value of the liquid phase fraction 5'^' as a 
function of the path parameter Xj"' 

VII. CONCLUSIONS 

In summary we have described a computational pro- 
cedure for accurately and directly determining the phase 
behaviour that occurs along a path in the space of the 
overall composition of a ternary fluid mixture. For each 
state point along the path the method yields estimates 
for the concentrations of the coexisting phases whose tie 
line intersects that point, together with the phase frac- 
tions. Since flnite-size corrections are exponentially small 
in the system size, the results are highly accurate even 



6 



when the path of interest touches the coexistence bound- 
ary where the fraction of one phase vanishes. The sim- 
ulation scheme, which is based on an extended samphng 
scheme in the semigrand canonical ensemble, uses his- 
togram reweighting to track the path of interest. The 
only requirement is to bootstrap the method by deter- 
mining the set of chemical potential differences that yield 
the overall concentrations corresponding to some point 
on the path. We have outlined how this can be achieved 
for pathway points corresponding to a single phase re- 
gion of the phase diagram, or within a coexistence region 
in the limit in which one component vanishes, ie. for a 
binary system. 



Our method is suitable for both fluid-fluid and fluid- 
solid phase coexistence, although in the latter case it 
needs to be combined with a sampling scheme such as 
phase switch Monte Carlo which allows the configuration 
spaces of the two phase to be connected by a trouble free 
sampling path [21, 22]. By virtue of the flexibility of his- 
togram extrapolation, it is also readily generalizable to 
paths that involve simultaneous variation in composition 
and an external field (temperature or pressure). In fu- 
ture work we plan to demonstrate its utility for tracking 
such paths in the contexts of mixtures of four or more 
components for which compact phase diagrams do not 
exist. 
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